% Set Up
framelist = dir("X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\filtered\*.jpg");
framelist = struct2table(framelist);
impath = 'X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\filtered';

sfPath = 'X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\contrast and sf\infant clutter code';
pyrPath = 'X:/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/orientation';

% framelist = dir("/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/filtered/*.jpg");
% framelist = struct2table(framelist);
% impath = '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/filtered/';
% 
% sfPath = '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/contrast and sf/infant clutter code';
% pyrPath = '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/orientation';
% low frequency cut-off, removes frequencies below this point in amplitude
% plot and fit
cutoffLF = 0.1; 
FreqPlot = 30;   % Number of freqs to plot in quick first view of amp spectrum (fig 3)
% cpi to cpd division constant
cpdCon = 41;  % vertical image size in degrees of Looxcie (480 pixels)

% get all files in selected folder, recursively
numImages = height(framelist);
 
% Pick one of these 3 window forms to apply to images
windowType = 'cos';
%% set up params
% read in first image to get dimensions of all of them
fn = fullfile(impath, framelist(1).name);
currImage = imread(fn);

% dimension first to make sure there is enough memory
imageDims = size(currImage);
imageDims = imageDims(1:2);

% struct to store parameters and data
params.impath = repmat(impath,numImages,1);
params.fileNames = repmat("",numImages,1);
params.numImages = numImages;
params.imageDims = imageDims;

% store fourier images and vectors
params.radialFreqMag = zeros(numImages,min(imageDims)/2);
params.rmsImageCon = zeros(numImages,1);
params.edges_full = zeros(numImages,1);
params.edges_center = zeros(numImages,1);
params.totalPower = zeros(numImages,1);
params.domOr = zeros(height(numImages),1);
%% run images through for loop
for fileCounter = 28101:numImages
    %fn = fullfile(impath,framelist(fileCounter).name);
    fn = fullfile(impath,framelist.name{fileCounter});
    currImage = imread(fn);
        
    imageDims = size(currImage);
    imageDims = imageDims(1:2);
    
    h = imageDims(1)/2;
    w = imageDims(2)/2;
    
    I2 = rgb2gray(double(currImage)./256);
    I3 = imcrop(I2,[w-200,h-200,400,400]);
    
    edgez1 = edge(I2,'canny',[0.10 0.25]);
    edgez2 = edge(I3,'canny',[0.10 0.25]); % need to crop edgez1 instead - diff edges here
    
    params.edges_full(fileCounter) = sum(edgez1, 'all');
    params.edges_center(fileCounter) = sum(edgez2, 'all');
    
    cd(sfPath);
    %addpath(genpath('X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\contrast and sf\infant clutter code'))
    addpath(genpath('/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/contrast and sf/infant clutter code'))
    
    if strcmp(windowType,'cos')
        windowWidth = round(min(imageDims)*0.8);
        imageWindow = cosine2D(min(imageDims),windowWidth);
        %need tukeywin - from signal processing toolbox
    elseif strcmp(windowType,'gauss')
        numGaussSd = 6;                         % number of sd of gaussian window
        gaussSd = min(imageDims)/numGaussSd;    % sd of window, in pixels
        imageWindow = dogauss2d(min(imageDims),gaussSd,0,gaussSd);
    else
        imageWindow = ones(min(imageDims));
    end

    % get the central square region of the gray image for fourier analysis
    Looxcie = mean(currImage,3);
    Looxcie = getcenter(Looxcie,min(imageDims));
    % window image (USE THIS ONE: Window in contrast rather than luminance)
    LooxcieMeanLum = mean(mean(Looxcie));
    Looxcie = ((Looxcie-LooxcieMeanLum).*imageWindow)+LooxcieMeanLum;
    % window image in luminance

    currFftImage = fft2(Looxcie);                   % take the fft
    currMagImage = fftshift(abs(currFftImage));
    currMagImage = currMagImage./min(imageDims)^2;      % divide by the number of pixels to scale it correctly

    % generate the summary amplitude spectrum
    ftmatrix =currMagImage;
    % make radii coordinates
    c = size(ftmatrix,1)/2+1;
    rmat = (1:size(ftmatrix,1))'*ones(1,size(ftmatrix,1))-c;
    cmat = rmat';
    radiusmat = round(sqrt(rmat.^2 + cmat.^2));
    freq = 0:max(size(radiusmat))/2-1;
    freqamp = zeros(size(freq));
    for i = 1:length(freq)
        currfreq = freq(i);
        currlocs = intersect(find(radiusmat>=currfreq),find(radiusmat<currfreq+1));
        freqamp(i) = sum(ftmatrix(currlocs(:)));  %% sum all of the energy at locations with that frequency
    end

    % output:
    params.freq = freq;
    % convert the frequencies from units of cycles per image (cpi) to cycles per degree (cpd)
    params.freq = params.freq./cpdCon;

    params.radialFreqMag(fileCounter,:) = freqamp;
    params.rmsImageCon(fileCounter) = sqrt(mean((Looxcie(:)-mean(Looxcie(:))).^2));
    % compute total power in fourier spectrum 
    params.totalPower(fileCounter) = sum(params.radialFreqMag(fileCounter,:));

    % calculate the slope of the function relating frequency to amplitude
    % (in log-log space)
    keeplocs = find(params.freq>=cutoffLF); % locations of frequencies to keep
    inputVals = [params.freq(keeplocs)',params.radialFreqMag(fileCounter,keeplocs)']; % x and y values for fit
    showFit = 0; % toggles plotting the data (1 = show, 0 = suppress)
    initialParams = [max(params.radialFreqMag(keeplocs)),1.5]; % initial parameter estimates
    [params.AmpSpecSlope(fileCounter,1:2), ~ ] =  invffitfn(inputVals,showFit,[],initialParams); 
    
    im = rgb2gray(currImage);
    nOri = 16; % number of orientations to perform the calculations on.

    % Some characteristics of the image.
    nPixels = numel(im);
    dims = size(im); % size of the image
    ctr = ceil((dims+0.5)/2); % center of the image

    % Calculate Fourier Transform and corresponding frequency axis.
    Fs = size(im,1);
    ff = -Fs/2:(Fs/2-1);
    imdft = abs(fftshift(fft2(im)))/nPixels; % amplitude (normalized by nSamples)

    % pre-calculate meshes for computing the filters
    [XX,YY] = meshgrid( ([1:dims(2)]-ctr(2))./(dims(2)/2), ...
        ([1:dims(1)]-ctr(1))./(dims(1)/2) );
    TH = atan2(YY,XX);
    rad = sqrt(XX.^2 + YY.^2);
    rad(ctr(1),ctr(2)) =  rad(ctr(1),ctr(2)-1);  % this line is getting rid of zero
    log_rad  = log2(rad);

    % pre-calculate position for the spatial frequency filter generation
    sf_positions = -[1/16,.5:.5:(round(log2(min(dims)))/2 +.5)]; % top/bottom edges of the filters
    sf_edges = size(im,1)*2.^(sf_positions-1);
    nSF = length(sf_positions)-1;

    % which orientations are we calculating this at.
    orientation = 90:-(180/nOri):-78.75; % Erin was right about the filter order and their orientations.

    % pre-allocate some storage variables for the loop.
    ssq = zeros(size(im));
    cnt = 0;
    amp = zeros(nSF,nOri);
    mid_freq = zeros(nSF,1);
    
    restoredefaultpath
    pyrPath = '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/orientation';
    cd(pyrPath);
    addpath('matlabPyrTools-master');

    % the loop
for ee = 1:nSF
    
        % RADIAL FILTER
        
        % upper bound
        if sf_positions(ee) == 1/16, width = 1/16;
        else, width = 1/4; 
        end
        position = sf_positions(ee);
        upper_radial_mask = generateRadialFilter(width,position,log_rad,true);

        % lower bound
        if ee < length(sf_positions)
            width=1/4;
            position = sf_positions(ee+1);
            lower_radial_mask = generateRadialFilter(width,position,log_rad,false);
        else
            lower_radial_mask = ones(size(upper_radial_mask)); 
        end
        
    radial_mask = upper_radial_mask.*lower_radial_mask;
    for aa=1:nOri

        % ANGLE FILTER
        anglemask = generateAngleFilter(aa,nOri,TH);
        mask = abs(anglemask).*radial_mask;

        cnt=cnt+1;
        
        % this can be used to demonstrate that the sum of squares of all
        % filters adds to 1 (except for at DC and high-frequency leftovers)
        ssq = ssq+mask.^2; 
        
        % selecting the amplitude associated with this filter.
        amp(nSF-ee+1,aa) = sum((mask(:).*imdft(:)));

    end
    
    % calculating the mid-point frequency associated with this SF band
    mid_freq(ee) = mean(sf_edges(ee:ee+1));
end

    amps_by_or = sum(amp);
    tot = sum(amps_by_or);
    amp_prop = amps_by_or/tot;
    A = sort(amp_prop,'descend');

    params.domOr(fileCounter) = A(1);
    
    if mod(fileCounter,100) == 0
               disp(fileCounter)
    end
end
%% output csv file
newDF = array2table(params.radialFreqMag);
cnames = strings(length(params.freq),1);
for i = 1:length(params.freq)
    cnames(i)=sprintf('%d',params.freq(i));
end

newDF.Properties.VariableNames = cnames;
newDF.fileNames = framelist.name;

writetable(newDF, 'X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\contrast and sf\sfvals_filtered2.csv')
%% SF - need to extract sf
low = NaN(length(params.domOr),1);
mid = NaN(length(params.domOr),1);
high = NaN(length(params.domOr),1);
for i = 28101:length(params.domOr)
   L = sum(params.radialFreqMag(i,1:80));
   M = sum(params.radialFreqMag(i,81:160));
   H = sum(params.radialFreqMag(i,161:240));
   tot = sum(params.radialFreqMag(i,1:240));
   low(i) = L/tot;
   mid(i) = M/tot;
   high(i) = H/tot; 
end

SF_sum = array2table([params.rmsImageCon(1:length(params.domOr)),params.edges_full(1:length(params.domOr))]);
SF_sum.Properties.VariableNames = {'rmsImageCon','edges_full'};
SF_sum.domOr = params.domOr;
SF_sum.filename = framelist.name(1:length(params.domOr));
SF_sum.low = low;
SF_sum.mid = mid;
SF_sum.high = high;

writetable(SF_sum, 'X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\contrast and sf\properties_filtered2.csv')

%% combined old and new
old = readtable('X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\contrast and sf\properties_filtered.csv');
for i = 1:height(old)
    if SF_sum.domOr(i) == 0
        SF_sum(i,:) = old(i,:);
    end
end
writetable(SF_sum, 'X:\1A_LooxcieBloomington\05_coding\13_VisualStatistics_2022\contrast and sf\properties_filtered_all.csv')